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Preliminary results of extensive numerical experiments with a family of simple 



["t i , models specified by the smooth canonical strongly chaotic 2D-map with global vir- 

\ > q \ tual invariant curves (VICs) are presented. We focus on the statistics of the diffusion 

rate D of individual trajectories for various fixed values of the model perturbation 
parameters K and d. Our previous conjecture on the fractal statistics determined 
by the critical structure of both the phase space and the motion is confirmed and 
studied in some detail. Particularly, we have found additional characteristics of 
what we termed earlier the VIC diffusion suppression which is related to a new very 
specific type of the critical structure. A surprising example of ergodic motion with 
a "hidden" critical structure strongly affecting the diffusion rate was also encoun- 
tered. At a weak perturbation [K <C 1) we discovered a very peculiar diffusion 
regime with the diffusion rate D = K 2 /3 as in the opposite limit of strong (K ^> 1) 
uncorrelated perturbation but, to the contrary, with strong correlations and for a 
very short time only. As yet, we have no definite explanation of such a controversial 
CN , behavior. 

o 

q ■ 1. Introduction: virtual invariant curves 

In a 2D-map (2.1) we are going to study here the diffusion crucially depends on 
the global invariant curves (GICs) which cut the 2D phase space (a cylinder, see next 
Section) of the motion. Even a single such curve is sufficient to completely block the 
global diffusion (GD) in the action variable along the cylinder. As is well known by now 
^ the existence of the GICs does depend not only on the perturbation strength but also 

on its smoothness. It is convenient to characterized the latter by the temporal Fourier 
spectrum of the perturbation. For analytical perturbation the Fourier amplitudes decay 
exponentially fast. In this case GD sets up if the perturbation e>e cr exceeds some critical 
value. Otherwise, the chaos remains localized within relatively narrow chaotic layers of 
nonlinear resonances. As a result, the GD is either completely blocked by GICs or the 
rate of the diffusion as well as the measure of its domain decay exponentially in parameter 
1/e as e — > (the so-called Arnold diffusion, for a general review see, e.g., |], 

By definition, the Hamiltonian of a smooth system has the power-law Fourier spec- 
trum with a certain exponent (3 + 1 (see, e.g., Q and references therein). In this case 
the GD is always blocked for some sufficiently small perturbation strength e < e cr (/3) 
provided the smoothness parameter (3 > (3 cr exceeds the critical value. This is similar 
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to the analytical Hamiltonian except that the critical perturbation depends now on the 
Hamiltonian smoothness (e cr (/3) — ► as (3 — > f3 cr ). 

To the best of our knowledge, the most strong rigorous result reads: f3 cr < 4 for a 
2D-map like in this paper (see ||). However, a simple physical consideration M leads 
to even smaller value (3 cr = 3 which is still to be confirmed somehow, theoretically or 
numerically. In any event, the smoothness of our model here /3 = 2 is even less. 

Until recently, the behavior of dynamical systems in the opposite case (3 < j3 cr of a 
poor smoothness remained rather vague. Even though the most numerical data seemed to 
confirm the simplest behavior of some universal GD (see, e.g,, ||) a few counterexamples 
were observed too (see, e.g., 0, ||). In the latter some trajectories remained within a 
certain restricted part of the phase space for a fairly long computation time. No clear 
explanation of these strange events has been given as yet. 

Meanwhile, about 20 years ago (!) a number of mathematical studies revealed var- 
ious possibilities for the existence of GICs in smooth systems with (3 < (3 cr (see, e.g., 
|10| , H |§). To us, the most comprehensive analysis of this problem was given by Bullett 
who rigorously proved a strange survival of infinitely many GICs amid a strong local 
chaos. Surprisingly, all these interesting results remain essentially unknown, at least to 
physicists. Apparently, this is because the abovementioned mathematical papers were re- 
stricted (perforce !) to what could be done rigorously that is to the invariant curves only 
without any attempt to analyse the very interesting and important transport processes 
such as diffusion. This is still in reach to the physical analysis and numerical (or labora- 
tory) experiments only. As a result, only after recent accidental rediscovery of GICs in 
chaos by Ovsyannikov [[Uj (which are still unpublished (!), see [12], [L3| for the full text 
of Ovsyannikov's theorem) the intensive physical studies of this interesting phenomenon 
have begun §0|I|§§. 

Interestingly, both authors |9|, [Tl| made use of exactly the same model for which a 
strange locked-in trajectory was observed still much earlier J/J. Apparently, this is because 
such a model (a particular case of our model with parameter d = 1/2, see next Section) 
is the simplest one possessing those curious GICs (for discussion see |T5|| ). Perhaps the 
main surprise was in that the GICs include the separatrices of nonlinear resonances which 
were always considered before as ones destroyed first by almost any perturbation. The 
principal difference is that now the invariant curves, separatrices including, do exist for 
the special values of the system parameters only (say K = K m ). 

Even though there are infinitely many such special values of the parameter, and the 
infinitely many GICs for each of the value of which just a single GIC does completely block 
the GD the probability of the latter that is the measure of such f^-values is apparently 
zero. So, a principal question to be answered is: what would be the behavior of that 
system for an arbitrary value of Kl In [03] we conjectured that even though the set of 
K m is not generally everywhere dense |J the density of this set is apparently rather high, 
so that one may expect some change (presumably suppression) of the diffusion for every 
K value as compared to the "usual" (familiar) dynamical system. In other words, we 
guessed that the structure of the phase space and of the motion therein may be changed 
by the formation of GIC at a close K value even if at almost any K there is actually no 
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GIC This is why we call now such a neighbour- K invariant curve the virtual one (VIC) 
with respect to any K |16| . 



Preliminary numerical experiments presented in [IB| did confirm our conjecture. This 
was done by the prompt computation of the average diffusion rate D(K) as a function of 
parameter K in the domain with GICs, real or virtual ones. The experiments revealed a 
very strong suppression of the diffusion, up to many orders of magnitude, restricted only 
by the computation time. But what turned out to be even more interesting was a very 
complicated (apparently fractal) structure of the dependence D(K). This seems to be a 
result of a very complicated structure of the model phase space itself. Preliminary, the 
latter looks like a so-called critical structure (see, e.g., H) but a rather specific one due 
to a forest of VICs. 

In the present paper we begin studying this seemingly new type of such a structure. 
Specifically, we start with investigation of the statistical properties of the diffuson as one 
of the characteristic processes in the chaotic motion. 

2. Model: the same again 

For reader's convenience we repeat below the discription of the model which is 
the same as in [15|, [16| . The model is specified by a map in canonical variables action 
(momentum) p - phase x: 

p = p + Kf(x) , x = x + p mod 1 . (2-1) 

where K = e > is perturbation strength (not necessarily weak), and "force" f(x) has a 
form of antisymmetric (/(— y) = —f(y), y = x — 1/2) piecewise linear "saw" of period 1. 
The phase space of the model is a cylinder: < x < 1 , — oo < p < +oo. 

As in fl"5| , 16] we are going to actually consider a family of maps with another parameter 



d (see Fig.l in JT5|) and the force 



/ 2x/(l -d), if \x\<(l-d)/2, 
nX) -\-2y/d, if|y|<d/2, ^ 

where y = x — 1/2, and the second parameter d (0 < d < 1) is the distance between the 
two "teeth" of the saw |/(x)| = 1 at points y = y± = ±d/2. The most studied particular 
case of the family corresponds to d = 1/2 when the two-teeth saw f(x) is symmetric. In 
the limit d = the two teeth merge in one, and all the invariant curves are destroyed. 
This was observed and explained in (15| for K > 0. In the opposite case K < (which 



is equivalent to K > , d = 1) the dynamics of the model is completely different, and 
we will not consider it in this paper (for a brief discussion see [|T^|). In our 2D-map (2.1) 



the GIC supports rotation of phase x around cylinder which bars any motion in p over 
GICs. Unlike this the local invariant curve (LIC), surrounding, for example, a domain of 
regular motion (see, e.g., and Section 5 below), corresponds to oscillation in phase x 
which allows other trajectories to bypass that obstacle. 

The GICs, separatrices including, do exist in the whole interval < d < 1 but for 
special K values only || 15, |16| . Particularly, the invariant curves are completely absent 
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for sufficiently large parameter 



2d 2 

K > K B (d) = , < d < 1 . (2.3) 

w 1 + d 

If K 3> i^B (see below) the physical quantity of the main interest for us, the diffusion 
rate D, can be approximately calculated from the Fourier expansion of force (2.2) (for 
detailes see |IJ) 

/(x) = £^sin(27mz), (2.4) 

n>l 71 

where 

_ 2 cos(nTr) sin(nTrrf) 
Particularly, in the limit d = 

2 

/n = cos(n7r) , (3=1. (2.6) 

7T 

the smoothness parameter (3 becomes less by one but the both values are less than critical 
f3 cr = 3. 

The calculation of the diffusion rate and other quantities are done using the standard 
analysis of the nonlinear resonances and their interaction (overlap) (see, e.g., [JX], |I5]|). 
The calculation is especially simple if one can neglect the variation of coefficients |/ n | ~ 
const in (2.4). This simplification is exact for d = (2.6), and remains reasonably accurate 



16] for 

K,3K B = ^ (2.7) 
Then the diffusion rate is approximately given by a very simple standard relation 



D(K) = ~ ^^5/2 w 57 ^5/2 (2 8) 
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where t is the motion time in map's iterations, and parameter K 1 is assumed to be 
sufficienly small. The latter expression in (2.8), which we will use below, is the result of 
extensive numerical experiments in confirmed also in [|l(| for K<0.1 (see |16[] and next 
Section). 

Notice that the dependence D(K) oc K 5 ^ 2 is different from the usual, or better to say, 
the simplest one D(K) oc K 2 . This is explained by the dynamical correlation of motion 
which is determined by the frequency of the phase oscillation on nonlinear resonances 



n 



n W a 



« 2^K « A n (K) « i (2.9) 



Here A n stands for the Lyapunov exponent charaterizing the local exponential instability 
of the motion which is the main criterion for dynamical chaos. Notice that for (3 = 1 
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both Q n and A n do not depend on the Fourier harmonic number n. The exact value of 
the Lyapunov exponent in the limit d = is given by 

A = In (1 + K + V2K + K 2 ) » V2K < 1 (2.10) 

The latter expression is the approximation for small K (cf. Eq.(2.9)) which is fairly well 
within the region of application of Eq.(2.8) (K<0.1) with accuracy ~ 1%. As the time 
in our model is discrete (the number of map's iterations) both correlation characteristics, 
(2.9) and (2.10), must be small, hence the above restriction on parameter K too. 

In the opposite limit K ^> 1 the correlation between successive x values are negligible, 
and one arrives at the "usual" relation for the diffusion rate: 

D(K) =K 2 C f\x)dx = % (2.11) 

JO 6 

independent of the parameter d. In intermediate region (K ~ 1) the correlation causes 
the decaying oscillation (see ||) which is beyond the scope of the present paper. 



3. Diffusion without any invariant curves: averages and moments 

As was already mentioned above there are no invariant curves for d = 0. Moreover, the 
motion is ergodic that is of the simplest structure of the phase space (cf. Section 4 below). 
Therefore, this particular case is not of the main interest to us by itself. Nevertheless, it 
is a good introduction to our central problem considered in Section 6 below. A similar 
approach was taken in our previous paper |16| . 



First, we consider the time dependence of the diffusion rate D(K;t). The semicolon 
here instead of usual comma is intended to emphasize that this time dependence is not 
a real physical contribution to the diffusion but rather a combination of two different 
processes: the proper diffusion via accumulation of random perturbation effects and a 
stationary regular oscillation of the diffusing variable (p in our case) which is a sort of 
background for the diffusion. Such a phenomenon can be roughly represented by a simple 
relation: 

D(K;t) ~ D^K) + (3.1) 
where B(K) is some function of the perturbation (see, e.g., [TBI and Eq.(3.5) below). In 



other words, in many cases, the present studies including, the nondiffusing stationary 
part can be separated from the diffusing part which much simplifies the analysis of this 
complicated process. All this can be describe, of course, via the standard method of the 
correlation of perturbation. However, this would lead to a much more intricate theoretical 
relations and, besides, to much less information on the diffusion dynamics (see, e.g., ||). 

An example of the diffusion kinetics is presented in Fig.l. The computation was done 
as follows. A number of trajectories M 3> 1 with random initial conditions homogeneously 
distributed within the unit area of the phase cylinder (0 < xq < 1 , < po < 1) were run 
for a sufficiently long time with successive outputs at certain intermediate moments of 
time t as shown in Fig.l. Remember that t is measured in the number of map's iterations. 
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Each output includes the diffusion rate < D > averaged over all M trajectories, and the 
dimensionless variance of that 



V < D 2 > — < D > 2 

For Gaussian distribution of action p this variance would be unity. This is indeed the case 
for a sufficiently long motion time when the measured diffusion rate reaches its asymptotic 
value -Dqo (3.1). A quite different dependence Vuit) for previous smaller time is of no 
surprise (nor is it very interesting) as over there D(t) depends on a completely different 
physical process one needs to pass over. 

A real surprise was the very beginning of the diffusion, the plateau in Fig.l. This 
looks as a real diffusion unlike the following part of a stationary oscillation. Moreover, 
the diffusion rate on the plateau Dq = K 2 /3 is the maximal one (2.11) as for large K 1. 
Another interesting observation is the duration of this strange diffusion 

1 1 

t() ~ — ~ — ; (o.o) 

A ^/2K 



which is close to the inversed Lyapunov exponent, the rise time of the local exponential 
instability of the underlaying chaotic motion. The last but not least curious property is 
the fast increase of variance (3.2): 

V M (t) ~ |, 2 < t<t (3.4) 

as shown in Fig.l. This is qualitatively different from the behavior of the same diffusion 
rate for large K 1 when the variance Vm ~ 1 is usual. The dynamical mechanism 
of this strange transitional diffusion is not completely clear and requires further studies. 
Apparently, it is related somehow to the main correlation (2.9) on the dynamical scale 
(3.3). Even though the initial "diffusion" is relatively fast it goes on a short time only so 



that the relative change of the initial distribution of trajectories |Ap|/|Ap| ~ \J D /A ~ 
K 3 / 4 <C 1 is negligible for small K <C 1 unless the initial distribution | Ap\ <K 3 ^ 4 is very 
narrow. However, in the latter case the dependence D(t) is very sensitive to the form of 
the initial distribution in p as a few our preliminary numerical experiments reveal. The 
variance of D(i) is especially strong for small t ~ t in the region of that mysterious 
plateau but eventually decays as t — > oo when the diffusion approaches its limit value 
.Doo. Apparently, this is related to a complicated fine structure of the phase space and/or 
of the motion correlations. This interesting question certainly deserves further studies but 
in the present paper we consider the simplest, homogeneous, distribution of the trajectory 
initial conditions on the phase cylinder. 

In this particular case, a very simple and surprisingly accurate empirical relation for 
the diffusion time dependence has been found starting from a qualitative picture (3.1). It 
has the form: 

Dp + tD 

(1 + T7) ] 



D(t) « /< ° - ^ , r = cAt (3.5) 
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Here r is the dimensionless time with an empirical fitting parameter c ~ 1 which is very 
close to one. The second empirical parameter 7 « 4 is less definite but it affects the 
turn of the dependence D(t) at r ~ 1 only. This relaxation of the diffusion rate has two 
time scales: (i) the plateau r v \ = 1 or t p i = 1/cA l/y/2K 3> 1, and (ii) the relaxation 
t r = D / Doo ~ 1/ ^> 1 or ~ 1/K which is still much longer. Interestingly, the 
usual diffusion spreading of a very narrow initial p-distribution on the relaxation time 
scale 

\Ap\ 2 R = D^tR = D^Do/D^/cA = D /cA = \Ap\ 2 pl 

is exactly equal to the spreading on the plateau. Hence, the full relaxation spreading is 
twice as large which is also directly seen from the empirical relation (3.5): 

1 nR v ; cA (1 + T lfh cA 

and which is still much less than the unit p-period. 

In Fig.l the empirical relation (3.5) is presented and compared to the numerical data 
in the dimensionless variables r and D* = D/D^ where is the asymptotic ("true") 
diffusion rate (2.8). In these variables the curves with various K values are similar, and 
converge in the limit r — > 00. 

Another interesting scaling can be done as follows. Let us calculate the diffusion rate 
Doo(D(r)) = D th from Eq.(3.5), and plot the ratio of that to the true rate (2.8): 

E* » • — + ^ ~ °" » 1 (3.6) 

Then, within the accuracy of scaling (3.5) and of fluctuations, this ratio must be always 
close to unity. This is indeed the case except the plateau (t<t ) where the rate D{r) is 
almost independent of r (see Fig.l). 

The next important statistical property are the fluctuations of the diffusion rate. One 
characteristic of those is the dispersion of tracjectories which is characterized by the 
variance (3.2). If all the trajectories would be statistically independent the dispersion of 
the mean diffusion rate were 

&<D> \ 2 2Vm 
<D> ) ~ M - 1 {6 ' n 

By construction, the trajectories are independent indeed with respect to their initial con- 
ditions but not necessarily to the corresponding diffusion rate. To check this we repeated 
the computation of diffusion N times with new and independent initial conditions, and 
then calculated the second (new) dimensionless variance for the average diffusion rate: 

(« D > 2 >n A M - 1 /o oN 

Again, if Eq.(3.7) holds true the variance Vn should be close to one. 
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The time dependence of both variances, Vm(£) and Vat(£)> is shown in Fig.l. Remark- 
ably, their behavior is qualitatively different. The first variance Vm{P) depends on the 
distribution function of p in the ensemble of trajectories while the second variance V)v(t) 
is affected by the statistical dependence (or independence) among trajectories whatever 
their distribution function. The results of our numerical experiments presented in Fig.l 
clearly demonstrate that the distribution in p quickly deviates from the Gaussian one 
during the diffusion on plateau, and come back only in the limit t — > oo when the dif- 
fusion rate D — > approaches the asymptotic value without any nondiffusing part. 
Unlike this, the trajectories remain statistically independent during the whole process of 
the diffusion relaxation. We will come back to discussion of this interesting point in the 
conclusion to this paper (Section 7). 

Now we turn to the most informative statistical characteristic, the distribution function 
f(D) of the diffusion rate. 



4. Diffusion without any invariant curves: the distribution function 

In the main part of our paper (Section 6) we shall be primarily interested in the 
distribution tail D — ■> of very low diffusion rate. The shape of this tail is known to be 
an important characteristic of the critical structure of the motion (see, e.g., @]). First 
indications of such a structure in the presence of the virtual invariant curves have been 



observed in ||16|| . Here we continue these studies. 

Since the statistics of the far tail is always rather poor we make use, as in fl6[], of a 
special version of the integral distribution 

F(D) = 1° f{D>)dD> « L (4.1) 
Jo J 



the so-called "rank-ordering statistics of extreme events" (see, e.g., fPT|). To this end the 
following simple ordering of the D(j) values (events) of the diffusion rate is sufficient: 
D(j + 1) > D(j), j = 1, 2, J. Then the integral probability is approximately given by 
the ratio j/J as shown in Eq.(4.1). 

In computation we typically ran M trajectories by iV times (see Section 3), so that the 
maximal number of events reached J = M x N = 10 4 x 10 = 10 5 . To obtain the lowest 
possible D values and minimize, at the same time, a rather big output we ordered all the 
computed events but printed out much less of those J J in such a way to get some 
all J\ < Jq first (the smallest) Dj while the rest were printed out in a logarithmic scale. 
An example of such distribution is presented in Fig.2 for K = 0.001 in variables D* = 
Dj < D > and F(D*) = j/J where < D > is some average diffusion rate (see below). 
The upper distribution corresponds to a rather long motion time t = 10 4 3> 1/K when 
the mean diffusion rate is already very close to the limit D^. For the lower distribution 
t — 10 is very short and corresponds to the plateau. 

At least in the former case when the p distribution is Gaussian (see Section 3) the 
distribution 

f(D) = f^D'-'e-^ (4.2) 
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is the so-called Pearson T-distribution with the two moments 



< D > = - , (AD) 2 =<D 2 > - <D> 2 = 4t (4.3) 
a or 

which are the mean and variance, respectively. For Gaussian p-distribution the reduced 
variance (3.2) Vm = 1 whence 

AD ' 2 ± = 2 (4.4, 



<D>J A 

and A = 1/2 independent of a. If, moreover, we introduce the dimensionless diffusion rate 

D - D* = JL (4.5) 

with average < D* >= 1, we obtain from Eq.(4.3) a = A = 1/2, too. Then, the new 
distribution becomes 

f(D*) = 



and 



2,71 



F(D*) = j° f(D')dD' -> \i — D i |.0) 




where the latter expression gives the asymptotics D* — > we need. This asymptotics is in 
a very good agreement with empirical data in Fig. 2 even at D* rs 0.1 (!). For very small 
D* the accuracy of the agreement is limited by the fluctuations due to a few remaining 
points. The smallest D* = 8.3 x 10 -11 corresponds to the estimate D* min ~ 1/J 2 = 1CT 10 . 
Since the distribution f(D*) in (4.6) is also Gaussian of \J~D* the integral F(D*) 



admits a very simple approximation found in JTS] 



1 - cxp J- D *{ 2) , D* > 1/2 



The relative accuracy \AF/F\ < 0.05 of that approximation is better than 5% in the 
whole range of F. Actually, the accuracy is even much better except a narrow interval at 
D* ~ 1/2. 

Thus, the upper distribution in Fig.2, which describes the real diffusion at sufficiently 
long motion time, is well in agreement with the available theory. This is no longer the case 
for the lower distribution on the plateau. In itself, this is not a surprise as in the latter case, 
unlike the former one, the measured diffusion rate is mainly determined by nondiffusive 
processes. However, a very interesting feature of this nondiffusive distribution is in that 
the exponent of the power-law tail remains exactly the same as if the p-distribution were 
again a Gaussian one. The simplest explanation, quite plausible to us, is in that the far 
tail still represents such a distribution which is a part of the whole distribution according 
to our original picture expressed by estimate (3.1). One immediate inference would be 
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decrease of the tail probability if we still use the same variable D* = D/D^. This is 
indeed the case according to the data in Fig. 2 ! 

A more difficult problem is the quantitative estimate of the distribution shift for 
the motion time t<l/K when the ratio < D* >=< D(t) > / > 1. This shift can 
be characterized either via the probability decrease by Rp times for a fixed D* or via 
the increase of D* itself by Rd times for a fixed probability. Notice that on the tail 
Rd = R F due to the square-root dependence (4.6). Characteristic Rd seems to us more 
preferable since it describes the shift not only of the tail but also (qualitatively) of the 
whole distribution F(D*). 

After some playing with the data we have found the following empirical relation for 
the tail shift: 

Rd(D*) « D a , (4.7) 

where the new diffusion ratio 

D,(r) « + 1 (4.8) 

and the fitted exponent a = 0.45. 

The philosophy behind this relation is following. We start with our original picture 
of a combined diffusive/nondiffusive process (3.1) which is almost our final choice (4.8). 
However, at the beginning we seemed to improve the original relation by inclusion of our 
surprising discovery, the plateau. Specifically, we tried to make use of Eq.(3.5) which is in 
a good agreement with the empirical data for the dependence D(t) (see Fig.l). Also, we 
have found that it partly described the distribution F(D) too, except on that mysterious 
plateau ! Then, our final, so far, step was the change from (3.5) back to a version of (3.1) 
in the form (4.8). 

How strange it may seem this did work with a reasonable accuracy as the insert in Fig. 2 
demonstrates. The remaining question "why?" is still to be answered in farther studies. 
Actually, this is a general serious problem of the dynamical mechanism underlaying the 
plateau formation and statistics. 

Our empirical relation (4.7) can be represented in a different way. Namely, instead 
describing the actual distribution tail shifted with respect to the limiting asymptotics 
(4.6) we may introduce the scaled diffusion rate 



D 

whence 



D 

Rd 



D* - ^- (4.9) 

The result is shown in Fig. 2 as a beam of 10 scaled distributions which are scattered now 
around asymptotics (4.6). 

5. Diffusion amid virtual invariant curves: the Lyapunov exponents 

Above we considered a very particular and most simple limiting case of our model (2.2) 
with parameter d — 0. In this case the motion is ergodic |§ which greatly simplifies the 
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problem under consideration. Nevertheless, we obtained a number of new results which 
form the firm foundations for further studies. 

The most important new feature of the motion for d > is the so-called divided phase 
space of the system that is a mixture of both chaotic as well as regular components of the 
motion. This is a typical structure of a few-freedom dynamical system (see, e.g., [f|). 

First of all, we need to exclude the regular trajectories from further analysis of the 
diffusion statistics. The standard well known method to do this is simultaneous compu- 
tation for each trajectory of the so-called Lyapunov exponent A that is the rate of the 
local exponential instability of the motion (see, e.g., p], [3| and references therein). In 
a 2D canonical (hamiltonian) map like our model (2.2) there are two Lyapunov expo- 
nents whose sum is always zero: Ai + A2 = 0. For a chaotic trajectory one exponent, 
say, Ai = A + > is positive while another one A2 = A_ < is negative. As a result, 
according to the standard definition of the Lyapunov exponent in the limit t — > 00, any 
tangent vector (dx, dp) of the linearized motion approaches the eigenvector corresponding 
to A + > 0. 

A simple well known procedure for computing A + , we made use in the present work too, 
is the following. For each of M trajectories with random initial conditions Xo,Pq we chose 
the tangent vector (dx, dp) of random direction and unit modulus: dp 2 = dx 2 + dp 2 = 1. 
Then both maps, the main one and the second one linearized with respect to the main 
reference trajectory x(t, xo,po), p(t,xo,Po) were run simultaneously during some time t. 
Finally, the current A(t) was calculated from the standard relation: 

A(t) = <lnp(t)> (5.1) 

t 

where the brackets denoted the averaging over M trajectories. Unlike the formal math- 
ematical definition of A in the limit t — > 00, in numerical experiments the Lyapunov 
exponent A(t) is always time dependent, perforce. 

In Fig. 3 a few typical examples of the A distribution are depicted for the number 
of events in (4.1) J = M equal to that of trajectories while the number of printed- 
out points Jo = M' < M is less except the case d = 0. The simplest one is for ergodic 
motion (d = 0). It has a form of almost vertical step which derivative dF/dA ~ 10 4 is very 
narrow 5-function. Notice that regular chain of points along F axis has no special physical 
meaning but simply reflects a particular accepted type of the distribution F(Aj) = j/J 
(4.1) with integer j. The mean A depends only on K (see Eq.(2.10)) but not on the 
initial conditions. This example in Fig. 3 shows the empirical/theoretical ratio which is 
very close to unity as expected. 

Two other examples correspond to the same K = 0.45 and M = 10 4 but different 
motion time t = 10 4 and 10 5 iterations. Both distributions have the same step at the 
largest A which corresponds to diffusive components (not necessarily a single one) of the 
motion similar to the ergodic case. However, the most interesting part of the former is 
the rest of the distribution which represents a rich motion structure contrary to a dull 
one in the ergodic motion. 

The largest (but again not the most interesting) part of this structure is related to 
the distribution steep cut-off at small A. Comparison of the two distributions for dif- 
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ferent motion time t = 10 4 and 10 5 shows that the A values of the trajectories in this 
region decrease with increasing time, approximately as A ~ 1/t. This would mean that 
all these trajectories are regular (see Eq.(5.1)) because the tangent vector p does not 
grow. The relative number of such trajectories gives the total area of regular motion 
on system's phase cylinder. In a particular example under consideration it amounts to 
A reg = 3177/10000 « 0.318 (t = 10 5 ). Generally, that value depends on a particular 
cut-off border from above chosen (see arrow in Fig.3). This delicate experimental prob- 
lem is considerably mitigated by a lucky feature of A distribution in our model, namely, a 
relatively wide plateau of F(A) immediately above the cut-off with only a few trajectories 
on it. However, the statistical accuracy 

A A 

^pl w {M . Areg) -m (5.2) 

is typically much worse, and can be improved by increasing the number of trajectories 
(and the computation time) only. 

Another interesting feature of A distribution in our model is a characteristic "fork" 
shape of the cut-off. This is a result of negative A for many regular trajectories. Such a 
peculiar representation is obtained by ordering A(t) values with their signs but plotting out 
the moduli |A(i)| only. Thus, the lower prong of the fork corresponds to negative A(t) < 
while on the upper one A(t) > are positive. This is due to the complex-conjugate 
Lyapunov exponents which results in a strictly bounded oscillation of the tangent vector 
(dx, dp) in this case. However, the area (5.4) A± ~ 0.20 < A reg « 0.318 is noticeably less 
than the total regular domains A reg . The rest is filled with the trajectories which are also 
regular but linearly unstable. This means the linear in time growth of the tangent vector 
p(t) ~ t so that A(t) — > remains positive but is vanishing in the limit t — > oo. This is 
the so-called marginal local instability with both A± = equal zero (for discussion see 
[T9|]). A curious point is that this seemingly exceptional case becomes the typical one 
in a nonlinear oscillator system due to the dependence of oscillation frequencies on the 
trajectory initial conditions. In fact, the bounded p oscillation producing negative A(t) is 
the exceptional case instead. The origin of this peculiarity is in a piecewise linear force 
in our model (2.2). As a result, the motion in the main (and the biggest for large K) 
regular domain around fixed point x — 1/2, p = is plainly a harmonic oscillation with 
the frequency (K < d) 

Q = arccos (l - —) » 1.47 (5.3) 



d J 

which remains the same in the whole regular domain of area 

2-kK 9 / K\ . , 

A - = I 1 - 2d) » 020 < 5 ' 4) 

Here y± = x± — 0.5 = ±d/2 is the position of two singularities of the force (see Eq.(2.2) 
and below) which restrict the size of the regular domain surrounded by the limiting ellipse 
to which both lines of singularity y± = d/2 = 0.25 are tangent. This ellipse is determined 
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by the initial conditions 

po = 0, x = 0.5 + y± (l - — ) « 0.5 ± 0.185 (5.5) 

All the numerical values above correspond to K — 0.45 and d = 1/2. Within the ellipse the 
motion of tangent vector obeys the same equation as the main motion, the only difference 
being an arbitrary length p of the tangent vector (for details see || and references therein). 

Coming back to Fig. 3 notice the decrease of the measured area A± with increasing 
motion time. This is explained by penetration of trajectories into a very complicated 
critical structure at the chaos border surrounding each regular domain (for details see, e.g., 
0]). For the same reason, the direct measurement of the whole regular region A reg « 0.40 
by a single chaotic trajectory as long as 10 9 iterations gives a noticeably larger value as 
compared to A reg ~ 0.318 from 10 4 trajectories for 10 5 iterations each. 

In all the curiosity of the A(t) distribution in regular components of the motion our 
main interest in the present study is the intermediate region between the regular cut- 
off at most small A(£) — > and the chaotic step at maximal A independent of t. The 
distribution in this region also does not depend on motion time and characterizes the 
proper critical structure of the chaotic motion. In an example in Fig. 3 this structure is 
presented by a relatively small probability step AF « 0.06 at A « 0.03. In the next 
Section a few other examples will be given too. 

6. Diffusion amid virtual invariant curves: the Critical statistics 

In Fig. 4 we present three characteristic examples of the effect of the critical structure 
on the diffusion statistics. The dashed curve shows the "unperturbed" distribution F(D*) 
(4.1) of the normalized diffusion rate D* = D/D norm where the normalizing rate D norm 
to be chosen in each particular case (see below). The term unperturbed means here the 
ergodic case d — without any invariant curves and critical structure (Section 4, the 
problem of critical structure in this case is not as simple as it may seem, see below and 
Section 7). In this case the normalizing rate D norm = is the true asymptotic diffusion 
rate (4.5). 

Now we are interested in the effect of the critical structure which typically arises in 
a nonergodic motion with its barriers for the chaos, or chaos borders. The latter are 
particular, and very important, case of an invariant curve which is transformed into itself 
under system's dynamics. As was discussed already above (Section 1) there are several 
different types of invariant curves (ICs). 

One is well studied and rather familiar chaos border surrounding any domain with 
regular motion. In this paper we call it the local invariant curve (LIC) which does not 
block the global diffusion (GD) around such a domain. An important property of LIC 
is the robustness which means that a small change of the system, say, of a parameter K 
or d cannot distroy the LIC but may only deform it slightly. This implies that LICs are 
always present in any divided phase space. 

Here we are mainly interested in a different IC type, the global invariant curves (GICs). 
Each GIC cuts the whole phase-space cylinder (x mod 1) of our model, and thus com- 
pletely prevents GD in p. Such ICs are less known, especially the most surprising of 
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them, the separatrix of a nonlinear resonance. However, those GICss are not robust in 
the model under consideration (see ||), being destroyed by almost any arbitrarily small 
perturbation of the system, particularly by a change of even a single its parameter. In 
other words, such GICs do exist for the special values, say, K = K m only. Even though 
there are typically infinitely many such special values the probability to find a GIC in a 
randomly chosen system is zero. This is why we are interested in a more generic situation 
when our model has no GICs at all. Yet, the effect of those still persists in a certain 
domain around each K m ! For this reason we call such GICs the virtual invariant curves 
(VICs) in analogy with other virtual quantities in physics like, for example, virtual en- 
ergy levels in qantum mechanics. Notice that unlike a GIC the VIC is robust and, hence, 
generic. 

Both LICs and GICs produce the so-called critical structure of the motion (see, e.g., 
H) which is typically characterized by a power-law distribution of principal quantities. 
The corresponding exponents c n are called the critical exponents. Their values are shown 
in Fig.4 at the related distributions. Notice that the opposite is generally not true that 
is a particular power law does not necessarily indicates any critical structure. In our 
model this is just the case for the ergodic motion where the diffusion rate distribution 
is also characterized by an asymptotic (D — > 0) power law with exponent Co = 0.5 (see 
above and Section 7). However, an important difference between ergodic and nonergodic 
dynamics is in that all critical exponents in the latter case c n < cq are less than (generally 
noncritical) ergodic exponent Cq. This is the main physical result of our preliminary 
numerical experiments we can present and discuss already right now (see Fig.4). 

Let us start with the distribution for K = 0.45 (upper solid line) which is far in the 
region without VICs (the border of this region is at Ks(d = 1/2) = 1/3, see Eq.(2.3) 
above and fll6| ). However, the regular trajectories (A reg « 0.318) together with LICs and 
the related critical structure are present. As a result the distribution (with D norm = D^) 
is well deviates from the unperturbed one for ergodic motion with d = 0. This type of the 
critical structure in a relatively narrow layer around a LIC is well studied by now (see, 
e.g., [|J) including some deviation of a typical distribution from a pure power law. The 
latter would mean the exact scale invariance of the underlaying critical structure in both 
the system phase space as well as its motion time. 

The critical structure is described by the so-called renormalization group, or renorm- 
group for brevity. On the other hand, the motion equations for any dynamical system also 
form a certain (dynamical) group. Such a fundamental similarity allows one to interpret 
the critical structure as a certain dynamics which was called the renormdynamics [^. [| . 
In this picture the exact scale invariance with a pure power-law distribution corresponds 
to the simplest, periodic, renormdynamics even though the original dynamics may be the 
most complicated chaotic motion. The resolution of this apparent paradox is in that the 
latter complexity of the original dynamics is "transferred" to the dynamical infinitely di- 
mensional space of the renormdynamics leaving behind the most simple renormdynamics 
itself (sometimes !). 

The latter limit is most studied simply because it is the simplest one. However, 
the generic case is just opposite that is a typical renormchaos is also chaotic f2T| , 
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Particularly, this implies a certain chaotic oscillation of the characteristic distribution 
around some average power law. This is just the case for the upper distribution in Fig. 4 
under discussion. It is characterized by the average critical exponent c\ = 0.3 with 
fluctuations ~ (c[ — ci) = 0.1. Such an interpretation of the critical structure in question 
is known to be typical but not necessarily unique (see below). The really unique property 
of this critical structure is the infinite power law, with or without the fluctuations. The 
term "infinite" corresponds here to the range of a renormdynamical variable lnD — > — oo 
with unrestricted variation even though the diffusion rate itself D > is strictly bounded 
from below. 

This is no longer the case for a new type of critical structure which we have encountered 
in our problem and which is produced by VICs (=robust GICs) rather than by the robust 
LICs. As was already explained above the principal difference between the two is in that 
the VIC is not an invariant curve at all. In terms of renormdynamics it means that VIC can 
mimic a GIC for relatively large In D only. This is clearly seen in Fig. 4 in the upper part of 
the distribution with local critical exponent c<i = 0.09 and parameter K = 0.335 (points). 
Here we have taken D norm = 10~ 6 < « 2 x 10~ 5 much smaller than the true diffusion 
rate D^. This shifts the whole distribution to the right to avoid the overlapping with 
other distributions. This value is slightly above the border K B (l/2) = 1/3 (see Eq.(2.3)) 
where there is a lot of VICs without any GIC. As a result the range of characteristic 
critical exponent C2 is very short: AlnD* « 5 as compared with the total available range 
« 25. The rest of the distribution remains fairly close to the unperturbed one. This 
would mean the absence of the critical structute over there or, at least, its sharp change 
at ln.D<2. In case of the former interpretation the renorm-motion stops in the latter 
region. 

In turn, this would imply a " dissipative" rather than "Hamiltonian" renormdynamics. 
Notice that the main part of the distribution is close but not identical to the unperturbed 
one including a slight difference in the characteristic exponent. Does it means a certain 
very slow renorm-motion remains a very interesting open question yet to be studied. 
Interestingly, the larger critical exponent c' 2 = 0.45 is also close to the local critical 
exponent c[ = 0.4 in the region without VICs or GICs and was interpreted above as a 
random fluctuation in the renormchaos. Is it really true still remains unclear. 

Finally, the third distribution in Fig. 4 (lower solid line) is actually coincides with 
the unperturbed distribution (D norm « Doo) even though it corresponds to the region 
with many VICs and strong suppression of the diffusion (K = 0.3294, see Fig.3 in ||16|| ). 
A deviation for very small D* is due to a poor statistics at this end. Notice that the 
coincidence of both distributions is not only asympotic (F — ► 0) but complete, including 
the opposite limit F — > 1. And this is in spite of a rather large regular region A reg « 0.581. 
The origin of this peculiarity for a particular K value remains unclear. One possibility 
is that, for some reasons, the area of the critical structure at the chaos border around 
this regular domain is unusually small. Examples of such a peculiarity in different models 



are known (see |[22|| ). Actually, in the latter work the critical structure was found to 
be unusually large but hidden. In other words, the motion was ergodic but with strong 
correlations (cf. the unusual diffusion rate (2.8) for K <C 1 in ergodic system at d — 0). 
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Coming back to this case in Fig.4 we can conclude that our "unperturbed" power-law 
distribution with exponent Co = 1/2 (dashed line) may well represent a peculiar critical 
structure related to the strong hidden temporal correlations rather than to a purely spatial 
geometry of the phase space. If this is true, indeed, the correlation decay may well be not 



a power-law one at all as is the case in model p2| where such a hidden decay is purely 
exponential (see Fig.6 over there). 

At last, let us mention another peculiarity of the critical structure in question: all the 
critical exponents found so far are less, if only a little, than "unperturbed" or "hidden" 
one Co = 1/2. The physical meaning of this universal inequality is in that the critical 
structure under consideration does always increase the probability of very low diffusion 
rate D —>■ 0. The general mechanism of this is known (see, e.g., Q), and is explained by 
the "sticking" of a trajectory within a complicated critical structure which slows down 
the diffusion. Interestingly, that the sign of the sticking effect may be opposite when the 
sticking accelerates the diffusion up to the absolute maximum D{t) oc t of homogeneous 



diffusion rate |4[] 



To summarize, we see that our "simple" model in the present paper really reveals a 
great variety of critical structure still to be farther studied and understood. 

7. Conclusion: A hidden critical structure ? 

In this paper we present some preliminary results of the numerical experiments with 
a family of simple models specified by the smooth canonical 2D-map (2.1) with global 



virtual invariant curves (VICs). As in our previous paper [16] we make use here of the 
same strongly chaotic model, and focus again on the statistics of the diffusion rate D which 
proves to be of a very complicated (apparently fractal) type determined by the so-called 
critical structure of both the phase space and the motion (see, e.g., ||). In previous paper 
|16H we studied the statistics of the mean diffusion rate < D(K) > averaged over ensemble 
of trajectories with random initial conditions. Our main result there was observation of 
very big and irregular fluctuations of the dependence < D(K) >, and a long and very 
slowly decaying tail of < D > distribution for < D >^ 0. We termed the latter effect 
the VIC diffusion suppression. 

In the present paper we continue studying this interesting phenomenon in more de- 
tails. To this end we turn from the statistics of averages < D(K) > as a function of 
model parameter K to that of individual trajectories for a given K. In principle, such 
an approach provides the most deep insight into a statistical problem. As the main sta- 
tistical characteristic we have chosen the integral distribution F(D) in the form (4.1) for 
a poor statistics at D — > 0. Preliminary results of our extensive numerical experiments 
presented in Fig.4 confirm, indeed, our earlier conjecture on a critical structure under- 



laying the fractal dependence < D(K) > in |L6|, the true sign of such a structure being 



various power-law distributions found. Moreover, besides the familiar well known critical 
structure exemplified in Fig.4 by the case with parameter K = 0.45 we did observe many 
cases of a rather different structure as one with K = 0.335. The principal difference of 
the latter is its finite size in the structure variable AlnD<5. The natural explanation of 
this difference is the following. First of all, the VIC is not a true invariant curve like GIC. 
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The latter completely blocks the GD while the former can, at most, inhibit the diffusion 
only. The mechanism of inhibition is known to be the sticking of a trajectory inside a very 
complicated critical structure. In turn, the sticking is the stronger (longer) the smaller 
is the spatial and/or the longer is the temporal scale of the critical structure. But for 
the VIC structure both are strictly restricted. On the other hand, such a restriction is 
the weaker the higher is the VIC density. In the system under consideration the latter 
is rather large, and hence the restriction leaves enough freedom for a strong suppression 
of GD for almost any K. Moreover, since the critical exponent of the VIC structure is 
typically very small (for example = 0.09 in Fig. 4) the probability of large suppression 
is high even for a short critical structure (cf. [16| for a different characteristic of this phe- 
nomenon). This slow decaying suppression probability is well ascertained in our numerics 
but we have not, as yet, any theoretical explanation of such a behavior. 

Now we come, perhaps, to the most interesting result of our current studies. How 
strange it may seem, this brings us to the apparently simplest case of our model with 
d = when the motion is ergodic. Can it still reveal any structure on the grounds that 
the distribution F(D) is also a power law (Fig.4)? That is the question! Certainly, it 
is not the case if, in addition, parameter K ^> 1 is big, and the diffusion rate has the 
standard form (2.11): D (x. K 2 . But what if K <C 1 is small? At least, the diffusion rate 
becomes qualitatively different: D cx K 5 ^ 2 . What does it mean? Generally, nothing! But 
in a particular case under consideration such dependence D(K) can be, and actually was, 
derived |16] from the resonance structure of the motion. If the system were not ergodic 
(with divided phase space) this structure would be clearly seen in the phase space. The 



question is what happens for ergodic motion with the same dependence D(K)? In [W 
we conjectured that some structure would still persist in the form of the correlations 
which determine the diffusion rate that is in some "hidden" form unseen directly in 
the picture of the motion in phase space. An example of such hidden critical structure 
was found in |E3 (see Section 6 above). However, in that case a particular distribution 



function was exponential rather than a power-law one (?). Hence, the question is if this 
qualitative difference could depend on a particular characteristic of the critical structure? 
Still another question arises from a very strange temporal behavior of the diffusion rate 
in the same "simple" case of ergodic motion for d = 0. We mean a "mysterious" plateau 
at the very beginning of diffusion under a weak perturbation (K <C 1, see Fig.l). In 
this case the dependence D(K) = K 2 /3 is the same as in the opposite limit of strong 
(K ^> 1) uncorrelated perturbation (?) but for a very short time only, the shorter the 
stronger is the perturbation (?!). Moreover, the correlations on the plateau are not only 
very large as in the weak-perturbation limit K — > but even increasing during the whole 
plateau regime (see Fig.l, dashed lines for variances Vm(t) (3.2)). As yet, we have no 
definite explanation for such a controversial behavior. A discreet current conjecture is 
the following. The duration of plateau is t p i ~ 1, or t p i ~ 1/A ~ 1/Q (see Eq.(2.9)). 
But the latter expression gives the period of phase oscillation on the critical nonlinear 
resonance which determines the diffusion rate fl6f . Then, one can image that this period 



does characterize not only the correlation decay, as usual, but also the correlation uprise. 
Still, the invariable diffusion rate over the whole plateau region is to be explained yet. 
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In any event, we are very curious to continue these most interesting studies! 
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Figure captions 

Fig.l Diffusion relaxation D*{r) = D(t)/Doo — > 1 in model (2.2) with parameter d = 
(without invariant curves) is presented as a function of dimensionless time r (3.5) 
for two values of K — 0.01 (circles) and 3 x 10~ 5 (crosses). Two smooth solid lines 
show the empirical relation (3.5) with two fitting parameters c = 1 and 7 = 4. 
Dashed lines are variances Vm(t) (3.2), and dotted lines show variances Vjv(t) (3.8) 
(see text). In the lower part the scaling (3.6) is presented reduced by factor 10 
to avoid overlapping with other data. The full volume of empirical data was J = 
M x N = 10 4 x 10 = 10 5 (see text). 

Fig. 2 Distribution function F(D*) (4.1) of the reduced diffusion rate D* (4.5) in model 
(2.2) without invariant curves (d = 0). Thick dashed straight line represents asymp- 
totics (4.6) of the integrated D-distribution (4.2) for the Gaussian p-statistics. Two 
lower wiggly lines correspond to large deviations from the latter: = 42 (K = 
10~ 3 ) and 461 (K = 3 x 10~ 5 ) (see Insert). A group of 10 D-distributions in a large 
interval (10 < < 461) are brought together using empirical relation (4.9). Insert: 
shift factor R D (see text) vs. deviation (4.8) for K = 10~ 3 (circles) and 3 x 10~ 5 
(crosses); the straight line is empirical relation (4.7). 

Fig. 3 Examples of distribution function F(A) of type (4.1) but for the Lyapunov exponent 
in model (2.2) for d = 0, M = M' = 80, t = 10 4 (the right-most step F(A), ergodic 
motion), and for d = 1/2, M — 10 4 , W = 1000, * = 10 4 , 10 5 (nonergodic motion, 
see text); in all cases K = 0.45. Horizontal line indicates the total share A reg m 0.318 
of the motion regular components. Arrow at A = 10~ 4 shows the lower border of 
chaotic trajectories chosen for further analysis (for t = 10 5 , see text). 

Fig. 4 Three characteristic examples of the diffusion statistics in the critical structure 
including virtual invariant curves (d = 1/2). Shown are integral distributions F(D*) 
(4.1) of the normalized diffusion rate D* = D/D norm (see text). The numbers at 
curves are critical diffusion exponents c m . The biggest one cq — 1/2 corresponds 
to the ergodic motion (d = 0) without any critical structure (dashed curve). Two 
straight lines show the average (c± = 0.3) and local (c[ = 0.4) critical exponents 
for K = 0.45 (solid line connecting 500 values of F(D*)). The distribution for 
K = 0.335 with two local critical exponents (c 2 = 0.09 and c' 2 = 0.45) is presented 
by 300 points shifted to the right to avoid overlapping with two other distributions. 
The third distribution (a solid line through 1000 points, K = 0.3294) is surprisingly 
close to that in ergodic case (dashed line). In all three examples M = 10 4 , t = 10 5 . 
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